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According to various embodiments, a method is provided for 
determining aberration data for an optical system. The 
method comprises collecting a data signal, and generating a 
pre-transformation algorithm. The data is pre-transformedby 
multiplying the data with the pre-transformation algorithm. A 
discrete Fourier transform of the pre-transformed data is per- 
formed in an iterative loop. The method further comprises 
back-transforming the data to generate aberration data. 

18 Claims, 2 Drawing Sheets 






U.S. Patent 


Oct. 9, 2012 


Sheet 1 of 2 


US 8,285,401 B2 



FIG.1 













US 8,285,401 B2 


1 

DISCRETE FOURIER TRANSFORM (DFT) 
ANALYSIS FOR APPLICATIONS USING 
ITERATIVE TRANSFORM METHODS 

STATEMENT OF GOVERNMENT INTEREST 5 

The invention described herein was made by an employee 
of the United States Government and may be manufactured 
and used by or for the Government of the United States of 
America for governmental purposes without the payment of 10 
any royalties thereon or therefor. 

FIELD 

The present invention relates generally to digital signal 15 
processing, and more particularly to a faster method of per- 
forming the discrete Fourier transform and the fast Fourier 
transform, in applications using iterative transform methods. 

BACKGROUND 20 

Fourier analysis is a family of mathematical techniques 
based on decomposing signals into sinusoids. The discrete 
Fourier transform (“DFT”) is the family member that can be 
used with digital signals and digital signal processing. A DFT 25 
can decompose a sequence or set of data values into compo- 
nents of different frequencies. In calculating each of the com- 
ponents (real and imaginary pairs) of the Fourier transform, 
one must use all N vectors of the sampled set. As is known in 
the art, the total computation time of the DFT is the total 30 
number of operations NxN, or N 2 . While the DFT is useful in 
many fields, computing it directly from the definition can 
often be too slow to be practical. 

A fast Fourier transform (“FFT”) is a way to compute the 
same result as the DFT, more quickly. Through the use of 35 
interlace decomposition, an FFT can compute approximately 
the same result as the DFT in N log N operations. Still this 
number of computations can often be slow and take up large 
amounts of processing time. A need exists to overcome the 
deficiencies of the DFT and the FFT. 40 

SUMMARY 

According to various embodiments, the present teachings 
provide an improved method for performing the discrete Fou- 45 
rier transform. By pre-transforming the data, performing the 
DFT, and then back transforming the data, the method 
described herein can compute the DFT in a total of N opera- 
tions. The present teachings describe the matrix form of the 
discrete Fourier transform (“DFT”) and the hybrid-index 50 
notation of Schouten is introduced as a method for perform- 
ing the DFT calculations in a complex vector space. The 
similarity transformation is presented in detail for a multi- 
dimensional DFT example (4x4 array size) and then applied 
to illustrate how the DFT can be carried out in a diagonal basis 55 
using the hybrid-index notation of Schouten. 

According to various embodiments, an optical system is 
provided that is adapted to collect optical data, for example, 
imaging data. The optical system can comprise a signal pro- 
cessor configured to carry out the DFT calculations described 60 
herein. By “pre-transforming” the data once at the beginning 
of the calculation, and “back-transfonning” the data at the 
end of the calculation, the DFT can be implemented as a 
one-dimensional vector product using just N operations, 
compared to the N log N operations needed by the fast Fourier 65 
transform (“FFT”) counterpart. The “pre-” and “back-” trans- 
forms are implemented using the DFT similarity transform 
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and the inverse of the DFT similarity transform. In some 
embodiments, the calculation can be implemented using arbi- 
trary frequency spacings. 

According to various embodiments of the present teach- 
ings, a method is provided for determining aberration data for 
an optical system. An optical system can be configured to 
perform the method. The optical system can comprise a signal 
processor, for example, one or more computers. The optical 
system can further comprise a control unit that can be oper- 
ably linked to the signal processor. The optical system can be 
launched into space prior to collecting data, for example, 
prior to collecting intensity data. The optical system can 
comprise a telescope. In some embodiments, the telescope 
can comprises one or more optical components, and at least 
one of the optical components can be configured to move. The 
optical components can comprise at least one mirror, lens, 
servo mechanism, prism, beam splitter, photomultiplier tube, 
detector, or a combination thereof. The control unit can be 
configured to control the movement of at least one of the 
optical components. The method can comprise collecting 
data generated by the optical system, for example, intensity 
data, light imaging data, thermal imaging data, or another 
type of data. 

According to various embodiments, the method can com- 
prise performing various calculations on the data. For 
example, a similarity transform, and an inverse of the simi- 
larity transform, can each be calculated based upon the data. 
In some embodiments, a discrete Fourier transform matrix 
can be calculated on the data. The method can comprise 
pre-transforming the data using the similarity transform and 
the inverse of the similarity transform to generate pre-trans- 
formed data. The method can perform an iterative transform 
operation on the pre-transformed data using the processor, by 
multiplying the pre-transformed data by the discrete Fourier 
transform matrix, for example, for a predetermined number 
of times, to generate iteratively transformed data. The method 
can back-transform the data using the signal processor, by 
multiplying the data by either the similarity transform or the 
inverse of the similarity transform, to form aberration data. 
The aberration data can be outputted to a user, for example, on 
a display such as a graph, and/or the aberration data can be 
printed out. 

In some embodiments, the method can comprise generat- 
ing correction data based upon the aberration data. In some 
embodiments, the optical system can be adjusted using the 
correction data, for example, an adjustment can be made to 
the orientation of an optical component under the control of 
the control unit. The aberration data can pertain to at least one 
of a piston aberration, a tip aberration, a tilt aberration, a 
defocus aberration, an astigmatism aberration, a coma aber- 
ration, a spherical aberration, and a trefoil aberration. 

In some embodiments, the iterative transform operation 
can comprise iteratively performing a phase retrieval opera- 
tion on the pre-transformed data. The method can further 
comprise generating phase correction data and the optical 
system can be adjusted using the phase correction data, for 
example, to overcome phase errors and provide an image that 
is in greater focus. For example, the correction data can 
comprise information that can adjust one or more optical 
components of the optical system to thereby compensate for 
aberrations in the data. 

According to various embodiments, the signal processor 
can be configured to generate correction data based upon the 
aberration data. The control unit can be configured to adjust 
the optical system based upon the correction data. The optical 
system can comprise a telescope, and the telescope can com- 
prise one or more optical components. The one or more opti- 



US 8,285,401 B2 


3 

cal components can comprise at least one of a mirror, lens, 
servo mechanism, prism, beam splitter, photomultiplier tube, 
detector, or a combination thereof. In some embodiments, the 
optical system can comprise a focusing system. In some 
embodiments, the optical system can comprise a beam-gen- 5 
erating device, such as a laser beam generating device, and a 
sighting system for the beam generating device. 

BRIEF DESCRIPTION OF THE DRAWINGS 

to 

The present teachings will be described with reference to 
the accompanying drawings, in which like elements are ref- 
erenced with like numbers. The drawings are intended to be 
exemplary only and not to limit the present teachings. 

FIG. lisa flow chart showing a DFT calculation method, 15 
according to various embodiments of the present teachings; 
and 

FIG. 2 is an illustration of an optical system comprising a 
telescope in communication with a signal processor, accord- 
ing to various embodiments of the present teachings. 20 

DETAILED DESCRIPTION OF THE PRESENT 
INVENTION 

According to various embodiments, the discrete Fourier 25 
transform calculation can be applied in diagonal form. The 
diagonal form introduces a transformation of basis by an 
application of the similarity transform of linear algebra. In 
doing so, a more efficient implementation of the DFT for 
applications using iterative transform methods can be 30 
achieved. The present teachings can be used in the field of 
signal processing, for example, the method can be used with 
phase retrieval for determining the phase error or aberrations 
in an optical system. 

Performance of optical systems, such as those found in 35 
cameras, telescopes, microscopes, and other systems for the 
collection of images, can often be improved by actively cor- 
recting for aberrations in the optical system. Wavefront sens- 
ing is a technique that can be used to obtain phase information 
from a wavefront to recover such aberrations. According to 40 
various embodiments of the present teachings, a method is 
provided for determining aberration data for an optical sys- 
tem. An optical system configured to perform the method is 
also provided. The optical system can comprise a signal pro- 
cessor. The optical system can comprise a telescope. The 45 
optical system can further comprise a control unit that can be 
operably linked to the signal processor, and, for example, an 
optical component of a telescope. The optical system can be 
launched into space prior to collecting intensity data. In some 
embodiments, the telescope can comprise one or more optical 50 
components, and at least one of the optical components can be 
configured to move. The optical components can comprise at 
least one mirror, lens, servo mechanism, prism, beam splitter, 
photomultiplier tube, detector, camera, or a combination 
thereof. The control unit can be configured to control move- 55 
ment of at least one of the optical components. The method 
can comprise collecting data generated by the optical system, 
for example, intensity data, light imaging data, thermal imag- 
ing data, or another type of data. 

According to various embodiments, the method can com- 60 
prise performing various calculations on collected data. For 
example, a similarity transform and an inverse of the similar- 
ity transform can be calculated based upon the data. In some 
embodiments, a discrete Fourier transform matrix can be 
calculated on the data. The method can comprise pre-trans- 65 
forming the data using the similarity transform and the 
inverse of the similarity transform to generate pre-trans- 
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formed data. The method can perform an iterative transform 
operation on the pre-transformed data using the processor, by 
multiplying the pre-transformed data by the discrete Fourier 
transform matrix, a predetermined amount of times, to gen- 
erate iteratively transformed data. The method can back- 
transform the data using the processor, by multiplying the 
data by one of the similarity transform and the inverse of the 
similarity transform, to form aberration data. The aberration 
data can be outputted to a user, for example, on a display such 
as a graphical display. 

In some embodiments, the method comprises generating 
correction data based upon the aberration data, and adjusting 
the optical system using the correction data. The aberration 
data can pertain to at least one of a piston aberration, a tip 
aberration, a tilt aberration, a defocus aberration, an astigma- 
tism aberration, a coma aberration, a spherical aberration, and 
a trefoil aberration. In some embodiments, the iterative trans- 
form operation can comprise iteratively performing a phase 
retrieval operation on the pre-transformed data. The method 
can further comprise generating phase correction data. The 
optical system can be adjusted using the phase correction 
data, to thereby overcome phase errors and provide an image 
that is in greater focus. For example, the correction data can 
comprise information that can adjust one or more optical 
components of the optical system to thereby overcome for 
aberrations in the data. 

According to various embodiments, the signal processor 
can be configured to generate correction data based upon the 
aberration data. The control unit can be configured to adjust 
the optical system based upon the correction data. The optical 
system can comprise a telescope, wherein the telescope can 
comprise one or more optical components. The one or more 
optical components can comprise at least one of a mirror, 
lens, servo mechanism, prism, beam splitter, photomultiplier 
tube, detector, camera, or a combination thereof. In some 
embodiments, the optical system can comprise a focusing 
system. In some embodiments, the optical system can com- 
prise a beam-generating device such as a laser beam-gener- 
ating device, and a sighting system for the beam-generating 
device. 

According to various embodiments, phase retrieval is just 
one example where the method of the present teachings can 
be applied but it will be appreciated that the method of the 
present invention is not limited to phase retrieval and can 
instead be utilized in various digital signal processing meth- 
ods . The method can be used anytime a transform of data from 
the time (or spatial) domain to the frequency domain is 
desired. The method of the present teachings can also be used 
anytime a transform of data from a frequency domain to the 
time (or spatial) domain is desired. In some embodiments, the 
method of the present teachings can be utilized in an optical 
system. The optical system can comprise one or more of a 
control unit, a signal processor, an optical unit, a telescope, 
and/or a combination thereof. In some embodiments, the 
optical system can comprise a space telescope, for example, 
the National Aeronautics and Space Administration’s 
(NASA’s) Hubble Telescope and NASA’ s James Webb Space 
Telescope (“JWST”). 

According to various embodiments, the optical system can 
be configured to collect optical data, for example, intensity 
data. The optical system can comprise an optical unit, for 
example, the fine guidance sensors of the Hubble Space Tele- 
scope. The optical unit can comprise one or more mirrors, 
lenses, servo mechanisms, prisms, beam splitters, photomul- 
tiplier tubes, detectors, cameras, or a combination thereof. In 
some embodiments, the optical unit can be disposed on a 
space telescope that is orbiting the earth. The optical system 
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and/or the optical unit can be manipulated by a user from the 
ground through the use of various wireless communication 
methods and wireless communication networks. A user can 
manipulate the optical system to adjust or move various opti- 
cal components within the optical system. This feature can 5 
allow for a user to adjust the optical system to overcome 
aberrations, obscurations, or other errors that can result from 
one or more of the components of the optical system becom- 
ing displaced, for example, during launch of a rocket carrying 
a space telescope that comprises the optical system. 10 

As is known in the art, to transform a digital signal from 
one domain to another, the discrete Fourier transform can be 
used. Diagonalizing the DFT matrix (the DFT eigen-prob- 
lem) can be a useful equation in transforming data in the field 
of signal processing. The matrix can be either one-dimen- 15 
sional or two-dimensional. Diagonalizing the DFT matrix is 
described, for example, in J. H. McClellan et al., “Eigenvalue 
and Eigenvector Decomposition of the Discrete Fourier 
Transform,” IEEE Trans, on Audio and Electroacoustics, Vol. 
AU-20 #1 (3), pp. 66-74, 1972, L. Auslander et ah, “Is com- 20 
puting with the finite Fourier transform pure or applied math- 
ematics?” Bull. Amer. Math. Soc., Vol. 1, pp. 847-891, 1979, 
and B. W. Dickinson et ah, “Eigenvectors and Functions of 
the Discrete Fourier Transform,” IEEE Trans. Acoust., 
Speech, Signal Processing, Vol. ASSP-30. #1 (2), pp. 25-31, 25 
1982, all of which are incorporated herein by reference in 
their entireties. 

Some aspects of the problem are related to topics in number 
theory that are solved by Gauss, as described in E. Landau, 
Vorlesungen uber Zahlentheorie , Vol. I. New York: Chelsea, 30 
1 927/1950 (reprint), p. 164, as described in L J. Good, “Ana- 
logues of Poisson’s Summation Formula,” Amer. Math. 
Monthly, Vol. 69, pp. 259-262 (1962), and J. H. McClellan, 
Comments on “Eigenvalue and eigenvector decomposition of 
the discrete Fourier transform,” IEEE Trans, on Audio and 35 
Electroacoustics, Vol. AU-21, p. 65, 1973, all of which are 
incorporated herein by reference in their entireties. For some 
applications, for example, phase retrieval, the diagonal struc- 
ture of the DFT is exploited in a special way, particularly 
when certain parts of the calculation do not have to be 40 
repeated at each iteration. 

Described herein are algorithms and equations that are 
useful in calculating a pre-transformation and a back-trans- 
formation of a digital signal. By pre-transforming the data, 
performing the DFT in diagonal form, and then back trans- 45 
forming the data, the discrete Fourier transform can be per- 
formed faster than the FFT, for certain applications. 

DFT as a Unitary Transform 

According to various embodiments, the DFT calculation 
can be implemented in matrix form, for example, as: 50 


N - 1 

F(v m ) s DFT(v m ) = /(*»)V™, 


where Y mn is a non-singular MxN Vandermonde array. 
Examples of the Vandermonde array are described in J. Tou 
“Determination of the inverse Vandermonde matrix.” IEEE 
Automatic Control (1964), p. 1820, and D. Balcan et al., 60 
“Alternative to the discrete Fourier transform.” IEEE Interna- 
tional Conference on Acoustics, Speech and Signal Process- 
ing (2008), pp. 3537-3540, both of which are incorporated 
herein by reference in their entireties. The Vandermonde 
array can be expressed as: 65 

y _ e -z2it v OT (n Ar) =e ~i2it(m n )IM = ( e ~2jtlM yn n =^ 1 J n 
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with components represented by: 
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Next, can be defined as: 

W mn = V mr /JM=(e- i2 * IM ) m n AfM, 

and then can be expressed as a unitary matrix: 

w-vft=i, 

using boldface matrix notation where I can be the identity 
matrix and the “ f” symbol can denote the complex-conjugate 
transpose or “adjoint.” The adjoint operation will be pre- 
sented more generally in the next Section below. Next, letting 

f„=Axf{x n \ 

the DFT calculation in matrix form takes the form of a unitary 
transformation: 

N - 1 

F(v m ) = >W». 

n = 0 

the inverse transformation can then be shown as (generally, 
taking M=N): 

M-l 

fn = Yj ( W ^ F m- 

m= 0 

Knowing the unitary matrix and the unitary transformation, 
the DFT can be analyzed using the transformation theory of 
complex vector spaces. The hybrid-index notation of 
Schouten supplies a natural framework for this approach and 
is described in J. A. Schouten, Ricci-Calculus, 2 nd edn., 
Springer- Verlag: Berlin (1954), and J. A. Schouten, Tensor 
Analysis for Physicists, 2 nd edn., Dover: New York (1959), 
both of which are incorporated herein in their entireties by 
reference. 

Hybrid Index Notation 

According to various embodiments, and as previously set 
forth, the DFT can be represented as a unitary transformation 
in a complex vector space. The calculations can be performed 
using the tensor-transformation formalism described by 
Schouten. Schouten’ s hybrid-index notation can provide a 
natural framework for the DFT calculations in a complex 
manifold. One goal in the kernel-index approach can be to 
give a more general representation of the DFT calculations by 
identifying components with quantities of a geometric mani- 
fold. Often, but not always, such identification can lead to 
alternative or more efficient calculation strategies that are 
discussed herein. By adopting a notation that makes the trans- 
formation rules of a quantity explicit, the meaning of the 
calculations can often be clarified and the notation can serve 
as an efficient bookkeeping device for checking consistency 
throughout the calculations. 
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DFT in Hybrid-Index Notation 

According to various embodiments, the DFT can be 
defined as a unitary transform between conjugate Fourier 
domains. The one-dimensional DFT transform can be 
expressed using Schouten’s hybrid-index notation as: 5 

F m =W m J n , 

and using the Einstein sum convention defined over a 
repeated indices in an upper/lower combination, the trans- 
form can be expressed as: 1C 




The upper/lower indices denote the covariant and contravari- 
ant transformation rules, respectively, for either a transforma- 
tion of coordinates, or a “point transformation” Allowable 
transformations of coordinates for this application will be 2 o 
discussed herein. The distinction between the transforma- 
tions of coordinates and point transformations are described 
in Schouten’s kernel-index method. To elaborate briefly, 
coordinate transformations are denoted by a transformation 
of basis, say, m— >m', where the new basis set is denoted by m\ 2 s 
The components in the transformation basis can be repre- 
sented as: 

u m '=A m m, u m , 

where A m 7m is the transformation matrix, the u m ' is the vector 30 
component in a new basis set, and u is arbitrary. A point 
transformation, meanwhile, can be represented as: 

U n =P n n u n , 

which is the form adopted for the DFT using the Schouten’s 35 
hybrid-index notation. Another distinction between a trans- 
formation of basis and a point transformation is that the 
indices for the point transformation are ordered and the posi- 
tions are marked using over- or under-dots as in P™„. The 
main take-away is that for coordinate transformations, the 40 
kernel symbol, u, remains the same. For point transforma- 
tions, however, the kernel symbol changes but the indices 
remain the same. The DFT transformation expressed using 
Schouten’s hybrid-index notation is defined as a point trans- 
formation, but when calculating the DFT in a diagonal basis, 45 
it is possible to consider the DFT using Schouten’s hybrid- 
index notation, such that one can use a common kernel sym- 
bol for the Fourier transform pair, for example, (F m ,I”)— > 
(f-,f 7 ). In some embodiments, the two transformations (basis 
and point) can be mixed, and can be used for transformations 50 
of basis in digital signal processing. 

According to various embodiments, the complex conjugate 
of the DFT using the hybrid- index notation can be expressed 
as: 

55 

where the over-bar can denote the complex conjugation of 
components and kernel symbols in accordance with 
Schouten’ s hybrid- index notation. Further, 

p>7i jyw 

The conjugation of indices and kernel objects can be used for 
calculations of the DFT transform. In some embodiments, the 
data values can be labeled by a contravariant component, f", 65 
and by complex conjugation of the kernel symbol and indices 
in various combinations, several other types of components 
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can be represented in Schouten’s notation by the raising and 
lowering of indices using the lundamental or metric 
tensor, a ki, for example, 

where the latter two can follow by complex conjugation of the 
first two. Contravariant component f" and its inverse can be 
defined as: 

and as mentioned above, covariant and contravariant trans- 
formation rules can be switched back and forth using the 
metric or fundamental tensor, a^, to raise or lower indices as 
appropriate. The metric tensor can define how the differential 
element of distance can be calculated in the geometric mani- 
fold as well as the inner product between two vectors, for 
example, u* and v l , thus, 

(u\v) =aj l u\ I =M / V / = 

The metric for raising and lowering indices can be repre- 
sented as: 


W™ = a ln W ml 

and 


r=^/ P 

and then substituting these expressions back into the DFT 
using Schouten’s hybrid-index notation, an alternative form 
of the DFT can be obtained: 

F m = w m n f n 
= 4 w "*fi 

=> F” = W^fn- 

Adjoint in a Geometric Manifold 

According to various embodiments, the adjoint operation 
applied to the W m „, denoted by **Adenoted by 


C = oo* 

can be defined through a two-step operation using both the 
metric tensor and complex conjugation, for example: 

(a) each index can be raised or lowered as appropriate using 
the metric tensor: 

(b) the result in (a) can be complex conjugated to give 

( ^J t= V E W k f W~™=w™, 

and therefore. 
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A comparison can be made in matrix boldface notation: 

(Wf=W=a Wa-^fT 1 , 

where “T” defines the transpose operation that is well known 
in the art. The matrix bold face notation illustrates a degree of 
non-generality in the way the adjoint can be calculated in 
matrix notation. Specifically, defining the adjoint using com- 
ponents in a geometric manifold (using steps (a) and (b) 
above) can give different results than complex conjugating 
and then performing the transpose operation using the “T”. In 
special coordinates, 

a nk = , 

these two can be identical, but in general, that is not always 
the case for arbitrary a- k . 

According to various embodiments, because the W m w , are 

unitary, the inverse, W m «can be defined using the adjoint and 
can be expressed as: 

k = k = ro + = wr, 

such that 

Considering the using two upper indices, W m>? , the cor- 
responding adjoint can be calculated: 

( W*y= a ^a ln W*=W^=W m -= W- m . 

The inverse DFT transformation to the hybrid-index notation 
can be defined as: 

W'nF™ = w n m w n 1 f=b n /=f, 

and when using Schouten’s hybrid-index notation, the DFT 
can be represented as: 

In some embodiments, in tensor transformation theory, 
there can be geometric meaning attributed to these differing 
but complementary transformation rules. Specific transfor- 
mations can be defined by an allowable group of transforma- 
tions, for example, G a . The corresponding allowable trans- 
formations of basis for the DFT application in digital signal 
processing, can include fractional and integer shifts of the 
data origin in addition to scaling and phase transformations. 
The details for these various transformations are discussed 
herein. 

Matrix Multiplication Using Index Notation 

According to various embodiments, when calculating with 
components, index order is important for proper matrix mul- 
tiplication and thus it is desirable to adopt the following 
convention for a rank-2 quantity: 

i 

col 


and the inverse of A m n , can be notated as, A 111 ^ such that 


where h m k is the Kronecker delta. To illustrate, the matrix 
multiplication of two matrices A and B can be notated as: 


and to be specific the matrix can be calculated using Math- 
ematica with 2x2 arrays. The Mathematic code can be repre- 
sented as, for example: 

Table[Sum [(A[[um, n]] B[[n, lk]]), {n, 1, NN}], {um, 1, 

2}, {lk, 1,2}]] 

where um denotes “upper-m” and lk denotes “lower-k,” the 
upper/lower dummy variable, n, can be left un-notated in 
Mathematica with regard to its position. The result can be 
represented as: 


a\b\ + a\b\ a\b; 2 + a; 2 tr 2 
a 2 { b\ + a 2 2 b\ a\b l 2 + a 2 2 b 2 2 


where AB is the conventional matrix product using matrix 
boldface notation. For comparison, reversing the order of the 
AB matrix to the BA matrix and transposing both the m-n and 
30 n-k index positions gives the BA matrix product: 


- A™ B? = BA = 


b\a\ + b, 2 cr^ b\a\ + b\ar 2 ' 
b\a\ + b 2 2 a 2 { b\ a\ 4 - b 2 2 a 2 2 ; 


and the corresponding Mathematica code for the BA matrix 
can be, for example: 

Table[Sum [(A[[n, um]] B[[lk, n]]), {n, 1, NN}], {um, 1, 
40 2}, {lk, 1,2}] 

Deriving the Metric Tensor from Unitary Invariance 

According to various embodiments, various geometric 
quantities can be defined that can be distinguished by their 
transformation rules using index position. For example, a 
45 group of allowable transformations, G a , can be introduced 
and then several invariant quantities can be constructed based 
on invariance of a given symmetry group. One such quantity 
can be the metric tensor, a ki, that can define distance in the 
manifold and can also define the inner product between any 
50 two vectors. The DFT metric tensor, a u, can be solved for 
postulating the unitary invariance under the action of the point 
transformation, W m „, as illustrated below. The can be 
defined as a unitary transformation and a metric tensor can be 
defined by invariance under the action of the W m n transfor- 
55 mation such that: 

or in matrix form: 

W a W~ i =a. 

60 Quantities that are invariant under the action of the unitary 
group can be hermitian and thus the a R can satisfy 


The introduction of a hermitian metric a B into the E„ geomet- 
65 ric manifold defines the manifold as a JJ„ which differs from 
the earlier notation used by Schouten, where this manifold 
was labeled as the FL. 
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General Form of ParsevaPs Theorem 

According to various embodiments, ParsevaFs theorem 
can be implied by unitary invariance, and not necessarily 
energy conservation as it is commonly referred to in the art. 
Consider the magnitude squared of F m — HFI 2 , that can be 
represented as: 

{F\F)=\F\ 2 =a-J 1 "F n . 

Substituting the DFT transformation rule from the hybrid- 
index notation equation we have: 


12 

The inverse transformation can be defined through the adjoint 

operation. Therefore, the inverse transformation, ,can be 
represented as: 

5 

W m n w\=8 m k , 

to and the =>can be given by: 


\F\ 2 =a- m W" k W{ff, 

but since the a B can be invariant under this unitary transfor- 
mation the right hand side can be simplified to: 


WT = w; 


_ 1 _ _J_ 

VaT V2~ 



\F\ 2 =a- k ff. 


Continuing to the general 2x2 metric, can be defined as 


ParsevaFs theorem is seen to hold for any solution of the a^, 
as amplified by unitary invariance, and thus in general: 20 

aJ fS F”=aJf. 


a It ~ 


' z Tl 

Z J2 ' 

( *11 + iyn 

*12 + iyn \ 


Z 12 , 

V *21 + *>21 

*22 + iyn ) 


The common form of ParsevaFs theorem can be expressed in 
vector notation as 

\F\ 2 =\f\ 2 . 


where the complex notation is dropped on indices of real 
components. The condition for unitary invariance then can 
provide the following simultaneous system of 4 equations and 
8 unknowns: 


In some embodiments, the standard DFT approach can be 
a special solution, for example, a coordinate system, such that 
for the a B given by the Kronecker delta: 


1 Z Jl + Z T2 + Z I1 + Z 12 Z Tl - Z T2 + %1 - %2 W z Tl Z T2 ' 

2 Z Jl + Z T2 ~ Z I1 ” Z 12 Z Tl ~ Z T2 ~ Z 21 + Z 22 Z 21 z 22 y 


a nk — &7ik i 

35 

and so in index notation ParsevaF s theorem in vector notation 
can take a Euclidean form: 


The unknowns of the system can have many possible solu- 
tions but will be further constrained by the heimitian condi- 
tion. Solving the system eliminates four of the unknowns (z Tl 
and Zj 2 ) in terms of the other four unknowns (z 21 and z 22 ) 
showing that a kl has the general form: 


or as it is commonly written: 

l-G i 2 +lz*2 i 2 + . . . +\F N \ 2 =\f l \ 2 +\f 2 \ 2 + . . . +\f N \ I 2 . 

According to various embodiments, the 5^ can be a solu- 
tion for the and can correspond to a positive definite sig- 
nature solution. The approach taken here can be motivated by 
that fact that in general, an eigen-problem requires the use of 
a general metric, and by no means does this necessarily cor- 
respond to the Kronecker delta in a unitary complex vector 
space, U„. Below is described how the Kronecker delta solu- 
tion is used to make progress on the main point of calculating 
the DFT in a diagonal basis, and for easy comparison with 
previously published results. 

Metric Solutions 
In a U 2 

By postulating unitary invariance of the metric one can solve 
for the The W m k can be identified with as The W m „, as 
previously set forth. Therefore, The W m k can be represented 
as: 


40 


45 


55 


f 2 %1 + Z 22 Z 21 i 
{ Z 11 z 22 ) 


Applying the heimitian (symmetry) condition to the other 
four unknown gives the condition that the 2x2 a Td can be 
valued: 


50 


( 2>21 + >22 

> 21 ) 

.f° °) 

V >21 

>22 ) 

do oj 


but there are multiple solutions for the a B in terms of the 2 final 
parameters x 21 and x 22 . The general solution for the 2x2 case 
can then be reduced from 8 to 2 unspecified real components 
and can be represented, for example, as: 


60 


_ ( 2*21 + *22 *21 ) 
V *21 *22 ) 





According to various embodiments, the solution can be 
checked by substituting the general solution back into the 
metric tensor. This verifies that this a kl , is invariant under the 
65 action of the for arbitrary, but real, numerical values of 

x 21 and x 22 . According to various embodiments, an infinite 
number of solutions can exist for the 2x2 metric, a H , as 
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implied by unitary invariance in terms of two real parameters. 
One possible solution is the special case of when x 21 =0 and 
x 22=1 giving the Kronecker delta, as mentioned earlier: 


14 

-continued 




The inverse DFT transformation can also be checked for 
consistency giving 


where the complex conjugation is kept for consistency, and 10 
the numeral “1” beneath the kernel symbol in Schouten’s 
notation labels this case as a distinct solution, different from 
other solution, say ? u. The labeling of kernel symbols in this 
way can be used in the eigen-problem discussed herein. 

According to various embodiments, and as mentioned pre- 
viously, there are an infinite number of solutions for the a^. 15 
Examples of such solutions can be the following, based on 
various choices for the x 21 and x 22 components: 





fl 

O'! 

1 


1 j 

x 2l^° ~ 

,0 

1 ) 

% \ *21 

= ll 

- 1 ) 

x 22^ 1 

x 22^~ l 





0 

1 


f 2 

n 

a kl 1 = 

~2 

l^l^ 1 = | 
x 22~*® 

lo 

1 } 

1 

1 



x 22^ 1 

U 2 


Using some specific solutions for the a^, one can examine 
various other quantities that are derived in Schouten’ s nota- 
tion. For example, consider the second solution above that is 30 
a multiple of the 


According to various embodiments, several other quanti- 
ties can be derived from the f” using the a B : 

f , =(z i ,z 2 );j^=a-j/ c (z l +z 2 ,z l -z 2 ) 

and comparing the right hand sides of the previous equations, 
as expected the result is f w =(f x ). Correspondingly, various 
other F m quantities can be given by 


F m ^-=(,z l +Z 2 ,Z l -Z 2 y, F m = aml F k = V2(z\z 2 ) 

VT 

r = J=(z T + Z 2 , Z T -Z 2 }. T m = ajm r = V 2 (z T , z 2 ). 


and components for the W mn and the W- w can be calculated: 




W mn = a kn W? = 


yfl 0 

0 V2" 


and then for generality define a complex data vector: 

f=(z l f)=(x l +iy\£+if). 

The DFT transformation can be represented as: 


In some embodiments, using the f-, and then comparing with 
40 the F m , one can check the forw ard DFT for consistency using 
the W mK : 


F m = w m f n = + ^ _ fy 

V2 


As a consistency check on both the solution for the a^, and the 
legitimacy of constructing the inverse transformation by use 
of the adjoint operation, i.e., 


VT ( Z l +Z 2 j 1,212 

W mn f n = . . , \=-=(Z l +Z 2 ,Z l ~Z 2 ) = F m . 

0 _i_ [z ~Z ) V2 

V2 


5 Q The adj oint operation on the W mn and W-„ in can be similarly 
verified: 


= (*0 = w n , 

one can check the inverse transformation calculated by way of 55 
the adjoint: 


j: _,)wC -Jm: .u 


<b& Wi a l7T 

i/i n 


V2U -1 


= w lk 

~ a mI a nlW m 


(i 1 W 2 if 1 1 

Hi -iJ 2I1 -1 

, V2 > 

( V2 0 1 

H 0 V2 [ 


such that 
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Finally, Parseval’s theorem can be demonstrated as: 

^^=aj//=lz‘l 2 -lz 2 l 2 , 

which is invariant. 

The General DFT 

According to various embodiments, further restrictions 
can be placed on the solution values for the x 21 and x 22 , by 


16 

According to various embodiments, the condition for uni- 
tary invariance can provide a simultaneous system of 16 
equations and 32 unknowns which can be further reduced to 
six free parameters by imposing the hermitian condition on 
5 the a%. After considerable algebra, the general solution for the 
a^, constrained by the hermitian condition, can be repre- 
sented as: 


a 


li ~ 


*33 + 2(x 4 i +* 43) 

*41 + 043 

(*33 + *41 “ *42 + *43 “ *44 + 2/}43) 
*41 + 043 


*41 - 043 
*44 

*43 - 043 
*42 


(*33 + *41 - *42 + *43 - *44 + 2l> 4 3) *41 - 043 \ 
*43 + 043 *42 

*33 *43 - 043 

*43 + 043 *44 


15 


classifying a s/ as being positive definite, negative definite, or 
indefinite, corresponding to the number of values (that are 
greater than zero, less than zero, or mixed, respectively), 
along the diagonal of the a#. When is definite there is no 
distinction between the various quantities that can be derived 
from &j d in the U 2 , and thus, numerical distinctions do not have 
to exist for components labeled by upper/lower indices, or 
even the complex conjugation of indices. Because a^ can be 
used to solve the eigen-problem, there can be different solu- 2 5 
tions for the eigen-problem corresponding to different solu- 
tions for the a**. Different solutions for the a^ can lead to 
different DFT eigen-quantities, but all the while preserving 
unitary invariance and the hermitian character of the a^. From 
the perspective of the transformation theory of complex vec- 30 
tor spaces, one can regard some properties of the DFT as 
deriving from the characteristics of the a#, for example, (a) a 
ki can be invariant under unitary transformations and (b) is 
hermitian. Condition (b) can be derived from (a) but it can be 
convenient when solving for the a e components, to regard (b) 35 
as an independent condition. 

In a U 4 

In some embodiments, the solution to a# corresponding to 
larger dimensional data sets can be solved for. Solutions 40 
where the are complex valued can provide consistency 
checks on the proper use of Schouten’s notation. By deriving 

~i 

a 4x4 W m k and its inverse, VF'%, that can be represented as: 


There are an infinite number of both real and complex 
valued 4x4 solutions for the a# in terms of the six free param- 
eters. Several example solutions corresponding to basic 
numerical values for the x and y components can be the 
following: 



( 1 

0 

0 

0> 


1 

0 

1 

0 

0 


1 * 33 -4, * 44 ->l 

0 

0 

1 

0 

’ 

*41->0,*42->0 

*43^°>3'43^° 

,0 

0 

0 

1 , 



f 

1 

—i 

2 i 

-/' 

1 


i 

1 

i 

0 

J \ ~ 

*33 -4, *44-4 
*41 -»0,*42-»0 

*43 ->0,3/43-4 


2 i 
i 

—i 

0 

1 

/ 

—i 
1 , 


*33^44 -* 1 
*41->-l+42->0 
*43 ->0,343 ->0 


1 -1 0 — 1 
-1110 
0 111' 
-1011, 


*33->l,*44~>l 

*41 ->-l,*42->0 
*43 O0/43 ->0 


-1 -1 -1 -r 
-1100 
- 1010 ' 
-1001, 


45 


'll 1 1 > 

1 1 -/ -1 / 
21-1 1 -1 
A / -1 -/, 


For exemplary purposes, a solution can illustrate a few of 
the basic quantities that are derived for the DFT in Schouten’ s 
notation. For generality, consider the following complex- 
50 valued solution for the a^: 



f 1 1 

1 1 ' 

1 

1 i 

-1 -i 

2 

1 -1 

1 -1 


A ~i 

-1 i , 


55 


1 -/ 2 / -r 
/ 1/0 
-2 i -i 1 -/ 

i 0 / 1 , 


and a general 4x4 metric can be defined as: and then ^ inverse can be represented as; 

60 


An 

^12 

^T3 

Z T4 ' 


( 1 -2-i 

-2 + 2/ -2 - / \ 

z ii 

Z 72 

Z 73 

Z 74 

-1 H * ^ 

a B =a =- 

-2 + 2/ 5 

2 + / 

-2 

%i 

%2 

%3 

%4 

1 

to 

1 

to 

to 

1 

1 

2-i 

An 

A2 

^43 

Z 44 y 

65 

(-2 + / -2 

2 + / 

5 , 
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To illustrate some specific quantities we can define a general ParsevaTs theorem can be invariant expected: 
complex “data” vector: 


f=C i r,Z,z 4 ), 


and then the DFT transformation can be represented as: 


a mn F m F n = a u f k f l = \z l \ 2 + \z 2 \ 2 + \z 3 \ 2 + \z 4 1 2 , 

As pointed our earlier, ParsevaTs theorem can take different 
forms for the different solutions. For example, consider: 


F m = W™f n = : 


z l +z 2 + z 3 + z 4 
z l - iz 2 - z 3 + iz 4 
z l -z 2 + z 3 -z 4 
z l + iz 2 -z 3 - iz 4 


Applying the inverse DFT transformation can show that: 


(11 1 1 \ 

' z l + z 2 + z 3 + z 4 N 

1 i -1 -i 

z l - iz 2 -z 3 + iz 4 

1-11-1 

A _2 . _3 _4 

z -z +z -z 

, 1 -i -1 i , 

1 , ■ 2 3 -4 


l Z +iz -z - iz J 


= (z l ,z 2 , z 3 , z 4 ) 

=/"■ 


a m nF m F n = a B f k f l = -\z l \ 2 + \z 2 \ 2 + \z 3 \ 2 + \z 4 \ 2 - 2 z l (z 2 + Z? + Z 4 ), 

which is invariant but not real valued. 

Example of Eigenvalue-Eigenvector Problem 
15 As previously set forth, multiple solutions for the a B exist 
based on unitary invariance, and thus, more general calcula- 
tions of the DFT and its various quantities can be extended to 
arbitrary basis sets. The additional generality stems from the 
possibility that the metric tensor does not necessarily have to 
20 be the Kronecker delta. As is known in the art, the eigenvalue- 
eigenvector problem (“eigen-problem”) seeks solutions to the 
following set equations: 


25 T mn v n = cr v m , 


Other quantities can be derived from the F” using the a#, for 
example: 


h ~ a Inf n ~ o 


Zl - izi + 2iz3 - iZA \ 
iZl +Z2 + iZ3 
-2 izi - iZ 2 +Z3- IZA 
izi + iz3 + ZA 


35 


and a further check can be performed to show that f x =(T w ), by 40 
calculating: 


Zl +iZ2~2iZ3 + IZA 


V 

for the vectors <;• of a general rank-2 quantity, T- r The scalar 

values, ?are the eigenvalues and the ). are the eigenvectors. 

v — 

The “a” index below the kernel object, in this case, a and t , 
does not refer to a component, but rather, the a-th 

eigenvector, \ , corresponding to the a-th eigenvalue, as 

separate and distinct solutions. Indices can be chosen from a 
set, for example, {a, b, c, d}, to label the kernel objects 
corresponding to distinct eigen-solutions, and then reserve an 
index range, for example, {k, 1, m, n}, to label the kernel 
objects corresponding to distinct eigen- solutions, and then 
reserve the index range, {k, 1, m, n}, to label components of a 
given basis. 

According to various embodiments, W- w can be substituted 
into the eigen-problem: 




-Hi +Z2-IZ3 
2 izi + iZ 2 + Z 3 + iza 
- iZl-iZ3+Z4 > 


45 


= 0 - v m , 

a a a 


Comparing f- to the complex conjugation of T w it can be seen 
that: 

/;= (a 

demonstrating consistency using the hybrid-index notation. 
Correspondingly, one can lower indices of the ¥ m to get the 
following result: 


thus switching the common kernel symbol for the Fourier 
transform pair, (F m , f 7 )— >(v-, v”), with the substitutions: 

and F^v-. 

The solution to the eigen-problem can result in a transforma- 
tion of the W- w into a diagonal basis by use of the similarity 
55 transform. Consider the overall structure of the DFT trans- 
formation in a diagonal basis for the 4x4 case: 



60 

WrrmV n = cr v m => 

■q 

o 

o 

o 



A/" 



r z l + (1 - 2 i)z 2 + (1 + 4 i)z 3 + (1 - 2i)z 4 ' 
(1 + 2i)z l - iz 2 - (1 - 2 i)iz 3 + iz 4 

1 

0 cr 0 0 
2 

' f ' 

f 


P 2 ' 

- cr 

r Fi'' 

f 2 > 

(1 - 4i)z l - (1 + 2 i)z 2 + z 3 - (1 + 2 i)z 4 

a a a 

o 

o 

o 

f 



a 

F V 

(1 + 2i)z l + iz 2 + (2/ - l)z 3 - iz 4 , 

65 

b^f 

o 

o 

o 

J 4 \ 
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where the primes on indices denote the data and its transform 
in the diagonal basis as discussed further below. The Fourier 
transform pair can be co-linear when calculated in the diago- 
nal basis: 


(t IzY 
{%lzY 


to 


(w mn - = 0 , 


20 

these degenerate eigenvalues are not necessarily orthogonal. 
The Gramm-Schmidt procedure as described in G. B. Arfken 
and H. J. Weber, “Mathematical Methods for Physicists,” 
Academic Press, 6 th . Ed., 2005, pp. 642-649, which is incor- 
porated herein in its entirety by reference, can be applied, or 
an orthogonal set can be obtained using the procedure 
described in R. Yarlagadda, “A note on the eigenvectors of 
DFT matrices,” IEEE Trans. Acoust., Speech, Signal Process- 
ing, Vol. AS SP-2 5, pp. 586-589, 1977, which is incorporated 
herein in its entirety by reference. 

Eigenvalues 

In some embodiments, for example, M=N=4, the DFT 
transform and its inverse can be represented as: 


that differs only by a scalar factor. This example illustrates the 
geometrical significance of the DFT transform in a diagonal 
basis and motivates using the same kernel symbol for the 
Fourier transform pair, since the quantities are co-linear. This 
difference in notation can be contrasted with the approach 
taken so far where the kernel symbols have been kept distinct, 
and the W m w transformation was regarded as a point transfor- 
mation. 

In some embodiments, the metric tensor can be used to 
show that the eigen-problem can be equivalently expressed in 
the form: 


15 


20 




n 

l 

1 

1 > 




l 

l 

—i 

-1 

i 



W.?: 

= 2 

l 

-1 

1 

-1 

; 




A 

i 

-1 

-i , 








f 1 

1 

1 1 

-1 




1 

1 

i 

-1 -i 

WT = 

= W 

T ~ 

:WT 









” 2 

1 ■ 

-1 

1 -1 






A 

-i 

-1 i , 


The solution can be solved for by: 


30 


or also as: 


det^W™ - <x<5™ j = 0, 


(w.™ - cr6 m n jv" = 0. 


35 


The eigen-problem does not have to depend on a specific 
choice for the metric. Solving the eigenvalues for the can 

be done by the following MATLAB® function: 


40 


after doing so, the 4x4 example has the solution of: 

I’cr- l^cr 4- l^cr + tj = 0 => jcr = 1, cr = -1, cr = 1, cr = -/J, 

with one degenerate eigenvalue, 


d et( VI 7 "' - cr<5!”) = 0. 

Eigen- Problem for the DFT Matrix 

According to various embodiments, it is known in the art 
that because 

W*=I, 

the eigenvalues for the DFT matrix for N>3 are non-unique 
and are multiples of the so-called “roots of unity”: 


<r = e tra/2 \a = i, ... , 4; => cr e {± 1, ±/}. 

a a 

Specifically, the eigenvalues become: 


cr = 1, cr = -1, cr = /, cr = -i. 

12 3 4 

The eigenvalues of the DFT transform are not all distinct for 
N^4, so there are many possible sets of eigenvectors satis- 
fying the eigen-problem. As a consequence of having these 
non-unique eigenvalues, the corresponding eigenvectors for 


cr — cr — 1 
1 3 

45 ... . 

Multiplicity of Eigenvalues for the DFT Matrix 
As discussed above, the significance of W 4 =I, for the DFT 
eigenvalue problem can be that degenerate eigenvalues 
appear for Ni=4. The multiplicities of these eigenvalues are 
5Q defined as the number of times each eigenvalue appears. For 
the 4x4 example, the result is 

m i =2, m 2 =l, m 3 =l , m 4 =0 , 

where the M* ordering corresponds to the ordering of the 
55 eigenvalues. These results can be generalized to higher 
dimensions as shown in the following Table of multiplicities 
for each eigenvalue: 


TABLE 1 


60 

Multiplicities for the DFT Matrix Eigenvalues 



mi 

m 2 

m 3 

m 4 


for cr = 1: 

for cr = -1: 

for cr = i: 

for cr = -i: 

M = N 

l 

2 

3 

4 

65 4m 

m + 1 

m 

m 

m - 1 

4m + 1 

m + 1 

m 

m 

m 
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TABLE 1 -continued 


22 


which can be simplified to the following: 


Multiplicities for the DFT Matrix Eigenvalues 


mi 

m 2 

m 3 

m 4 

M = N 

for cr = 1 : 
l 

for cr = - 1 : 
2 

for cr = i: 
3 

for cr = -i: 
4 

4m + 2 

m + 1 

m + 1 

m 

m 

4m + 3 

m + 1 

m + 1 

m + 1 

m 


to 


-v 1 + v 2 + v 3 + v 4 
1111 

v 1 - (2 + i)v 2 - v 3 iv 4 
l ill 

„1 ,,2 ,,3 ,,4 


v 1 + IV 2 - v 3 - (2 + i')v 4 
111 1 


The 4x4 matrix corresponds to the m=l value in the top row 
(4 m ) of the table giving the multiplicities of eigenvalues 15 
listed. Similarly, it can be shown that for the 5x5 case, the 
eigenvalues comprise the following solution: 


The solution for j can be given in terms of two arbitrary 
components, for example, V and / , thus leaving: 


20 


(v\v 2 , v 3 , v 4 ) = (v 3 +2v 4 , v 4 , v 3 , v 4 ). 
Vi l l l / Vi 1111 / 


j'cr- l^cr + l^cr 2 + 1) = 0 => jcr = 1, <x = -1, cr = /, cr = cr = 1 j, 

corresponding to the 2 nd row (4 m +l) of the Table with m=l, 
and M=N=4m+ 1 . Substituting in the value of 1 form, M=N=4 
(1)+1=5. Thus, the multiplicities for the 5x5 case are 

OT 1 =2 j OT2=1,/W3=1 j OT 4 =1. 

In another example, say M=N=131 : =(4x32)+3, correspond- 30 
ing to the fourth row (4 m +3) of the Table giving the following 
multiplicities: 


In this example, two components can be used because of the 
two-fold degeneracy of the 


cr = cr = 1 
1 3 


mi=33 J m ? =33M?=33,m 4 =32. 


35 


Eigenvectors 

According to various embodiments, given the eigenvalues, 
the eigenvectors corresponding to each eigenvalue can be 
constructed using the eigen-problem. It is instructive to con- 40 
sider this process for the eigenvectors corresponding to the 
degenerate first and third eigenvalues, 


cr = cr = l. 
l 3 


45 


eigenvalue. The eigenvectors corresponding to unique eigen- 
values can be specified in terms of a single component. Values 

3 4 

for V and I can be arbitrarily chosen, for example, 


v 4 = 1 and v 3 = 0, 

l l 


to obtain a realization for the first eigenvector: 

(v l ,v 2 , v 3 , v 4 Wv 3 +2v 4 , v 4 , v 3 ,v 4 )| 4 , 

VI 1 1 1 / VI 1111/ Ip* 1 1 

vVo 

1 

= (2, 1, 0, 1). 


The first eigenvector can be obtained as a solution to the set of 
linear equations: 


l=°= 

\ a 'a lQ~»l 


(l-o- 1 1 1 \ 

l 

v 1 

l 

1 -i-cr - 1 i 

l 

v 2 

1 

l 

l 

l 

v 3 

1 

, 1 ' - 1 - ; -T, 

v 4 
V 1 


50 A solution for 3 , can be obtained by choosing complemen- 

3 4 

tary but arbitrary numerical values for the and . , for 
example, 


55 


v 4 = 0 and v 3 = 1 : 
3 3 


60 


(v 1 , v 2 , v 3 , v 4 ) = (v 3 + 2v 4 , V 4 , v 3 , v 4 ) I 4 n => v" 

\3 3 3 3 / V3 3 3 3 3 / r 3 


65 


= ( 1 , 0 , 1 , 0 ). 
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In some embodiments, the remaining eigenvectors that 
correspond to their eigenvalues, 


cr = — 1 and cr = -i, 
2 4 


24 

-continued 
1 l-ll' 

-i , l-ll 1 1 

s ™ = - 

m 4 1-13-1 

,0 -2 0 2 , 


can be solved for. An exemplary set of eigenvectors (but not 
the only set) for the 4x4 DFT as transformation can be rep- 
resented as: 


It can then be verified that applying the similarity transfor- 
mation to the gives: 



f T 


f -1 ' 


f 1 ' 


( 0 ' 


1 


1 


0 


-1 

v n = 

0 

, v w = 


, v" = 


, v” = 

0 

1 

2 

1 

3 

1 

4 


,1; 


, 1 , 


,0, 


, 1 , 


According to various embodiments, because two of the 
eigenvalues are degenerate, neither the eigenvector v^, or the 
set of eigenvectors v 1 ' , -v 4 ”, form an orthogonal system that 
can be verified by forming the inner product between any two 
vectors using the base solution, 


a Ti = 6 Ti . 


where the off-diagonal entries correspond to the 30 


cr — cr — 1 
1 3 


35 

value. For reference, we can construct an orthonormal set 
from the set of eigenvectors: 


'2' 


r-r 


( 1 > 


( 0 ' 

1 

l 

l 

i 

-1 

1 

-1 






,«" = — 


0 

2 2 

l 

3 2yf3 

3 

4 sf2 

0 



, i , 


,-l , 


, 1 , 


and verify that 


5-rfltu 1 = 6 mn . 
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Similarity Transform 

Given the eigenvectors, one can form the similarity transform 
using the eigenvectors as columns. This can be denoted as a 
transformation of basis such that one can find by construction 

that the similarity transform, S m mt , and its inverse, S 2 can be 
given by: 


s™ = 


2-11 O' 
110-1 
0 110’ 
110 1, 
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, -l 

SWS~ l = S% W™S n n r 
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( cr 0 0 0 3 
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0 0 <r 0 

3 

0 0 0 <r 

4 ) 

n 0003 
0-100 
0 0 1 0 

,0 0 0 

As a check on index ordering one can run the following line of 
Mathematica code with Si defined as the inverse of the S: 
Exemplary code is listed below. 

Table [Sum[S[[m, ump]] W[[m, n]] Si[[lnp, n]], {n, 1, M}, 
{m, 1, M}], {ump, 1, M}, {lnp, 1, M}] 
where ump denote “upper-m-prime” and lnp is “lower-m- 
prime.” The similarity transform can be constructed using the 
orthonormal eigenvectors as columns. Where the similarity 
transform leaves the W m fl in diagonal form, many possible 
choices can be made for the S w mt . As indicated by calculations 
thus far, it does not seem to matter whether or not the eigen- 
vectors forming the S„ mt are orthogonal, and this can greatly 
simplify the initial calculations needed for calculating in the 
diagonal basis. 

DFT Calculation in Diagonal Form 

According to various embodiments, the diagonal form of 
the DFT can allow the DFT calculation to be performed as 
simply a multiplication of the entries on the main diagonal 
(the eigenvalues), to the corresponding data components that 
have been pre-multiplied using the similarity matrix. This 
means that after pre-multip lying the data, the DFT transform 
can be reduced to a vector product consisting of N operations 
rather than a matrix multiplied with N 2 operations when the 
Fourier transform is used, or N log N operations as is the case 
when a fast Fourier transform (“FFT”) is used. This can result 
in a reduction in computation time and cost. In some embodi- 
ments, the diagonal DFT representation can lead to a fast 
calculation for one part of the DFT operation, but this gain can 
be offset by the added computational cost of having to pre- 
transform the data and then back-transform to bring the Fou- 
rier transform into the original basis. In some embodiments, 
however, these latter operations need not be performed at 
every step in an iterative transform application, such as, for 
example, in a phase retrieval operation. 

According to various embodiments, once the notation of 
the DFT calculation in the diagonal basis has been switched 
back to distinct kernel symbols for the data and its transform, 
the DFT can be represented as: 


F m, =W m ' n f\ 
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where the m' represents the transformation of basis under the 
action of the similarity transform. The only non-zero values 
for the W m „ are along the diagonal, therefore, one can form a 

vector, 5 fro m the eigenvalues, for example, from the diago- 5 
nal entries of the W m n . Continuing with the calculation of the 
starting data, f k , can be pre-transfonned or pre-multiplied 


using the inverse similarity transform S„ 5 to give 


to 


/”' =~s n 'f n . 


The diagonal DFT can then be the product of the N compo- 15 
nents of the Sand the f„* defined, as shown: 



20 


where no sum is implied over the (m 1 , m') pair. To complete 
the calculation, the “diagonal” DFT components can be trans- 
formed back to the original DFT basis, F 7 ". To obtain this 25 
result one can express the left hand side of the DFT calcula- 
tion in the diagonal basis in terms of components in the 
original basis: 


p'tn S m 


Therefore, the original Fourier transform components can be 35 
obtained by back-transforming using the similarity trans- 
form, and be represented as: 


F n = s n n , F n> = s n n , s£ F m = 6 n m F m . 

The inverse DFT can be calculated in a diagonal basis because 
the inverse of a diagonal matrix can be found by populating its 
diagonal entries with the inverse of the individual compo- 
nents. Therefore, given the W m w , the inverse can be (no sum is 
implied on the right hand side over the (m', m') pair): 


</=( 

and then inverse DFT transformation in the diagonal basis can 
be represented as: 


60 

According to various embodiments, where it is feasible to 
perform the DFT in a diagonal basis, the similarity transform 
and its inverse can be pre-calculated. The multiplication of 
the similarity transform with the data can also be pre-calcu- 65 
lated. Then, at an appropriate point, the DFT results can be 
transformed back to the original basis using the inverse simi- 


larity transform, preferably at the end of an iterative loop. An 
exemplary algorithm approach can comprise the method 
illustrated in FIG. 1 . 

As shown in FIG. 1, in a first step, the similarity transform 
and the inverse similarity transform can be pre-calculated. In 
a second step, the diagonal form of the DFT matrix can be 
calculated in a third step wherein a data set can be pre- 
transformed by multiplying the data by the inverse similarity 
transform. In a fourth step, the DFT and the inverse DFT can 
be calculated for a set of data. This fourth step can be for an 
iterative transform application, for example, phase retrieval. 
The fourth step can loop a desired number of times, for 
example, a fixed number of times set by a user. 

FIG. 2 illustrates an example of an optical system 10 com- 
prising a telescope 20 and a control unit 22. Telescope 20 can 
comprise one or more optical components, for example, at 
least one mirror, lens, servo mechanism, prism, beam splitter, 
photomultiplier tube, detector, camera, or combination 
thereof. In the embodiment shown, telescope 20 can comprise 
a lens 34. Telescope 20 can be configured to be launched from 
the ground into the earth’s atmosphere. Telescope 20 can be 
configured to orbit the earth. Control unit 22 can comprise a 
signal processor 24. Signal processor 24 can be configured to 
carry out one or more of the algorithms of the present teach- 
ings. In some embodiments, one or more of the algorithms of 
the present teachings can be coded in software that can be 
processed by signal processor 24. Control unit 22 can be 
configured to control at least one of the one or more optical 
components of the telescope 20, for example, lens 34, through 
a wireless communication. In some embodiments, control 
unit 22 can also comprise a land-based transmitter/receiver 
26, and a space-based transmitter/receiver 28. Signals trans- 
mitted by transmitter/receiver 26 can be received by transmit- 
ter/receiver 28 and used by control mechanics 30 in telescope 
20, to control movement of lens 34. Control mechanics 30 can 
comprise appropriate motors, linkages, and gears to effect 
fine movement of lens 34. 

According to various embodiments, the algorithms 
described herein can be used with any desired signal process- 
ing application where a transfer of data between the spatial 
domain and the frequency domain is used, or vice versa. For 
example, the methods provided herein can be used with a 
hybrid diversity algorithm that is used to recover phase 
retrieval that has high spatial frequency. The methods pro- 
vided herein can be used with a distributed transpose algo- 
rithm for a two-dimensional FFT on a cluster of digital signal 
processors. In some embodiments, the methods provided 
herein can be used with matrix DFTs for fast phase retrieval. 
The methods herein can be used with space telescope wave- 
front sensing software, for example, the software used with 
NASA’s James Webb Space Telescope. In some embodi- 
ments, the methods provided herein can be used with hybrid 
diversity methods that utilize an adaptive diversity function as 
described in U.S. patent application Ser. No. 11/469,105, 
filed Aug. 31, 2006, to Dean, which is incorporated herein in 
its entirety by reference. The methods of the present teachings 
can also be incorporated with various digital signal process- 
ing applications for example, algorithms using the Nyquist- 
Shannon sampling theorem or variations thereof, algorithms 
that provide optimal padding for two-dimensional FFTs, and 
hybrid architecture active wavefront sensing and control. 

EXAMPLE 

The steps above are illustrated for a 4x4 example. Consider 
the general complex data vector: 

/„=(z\zV,z 4 ). 
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The inverse similarity transform times the data vector is given 
by: 


Zl +Z2-Z3+Z4 ) 


r' 



-Zl +Z 2 +Z 3 +Z 4 
Zl - Z2 + 3^3 - Z 4 


2(-z 2 +z 4 ) 


O’ = (1, -1, 1, -0- 


5 


10 


The DFT in the transformed basis is followed by multiplica- 
tion of the eigenvalues with the f"' composition to give the 
Fourier transform in the diagonal basis: 


Zl +Z2 ~Z3 +Z4 ) 


F m> — (T ■ f m> 


1 ZI-Z 2 -Z 3 -ZA 
4 Zl - Z2 + 3Z3 - ZA 
, -2i(-Z2+Z 4 ) 
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Finally, the f"* components can be back-transformed by mul- 25 
tiplying with the similarity transformation to give the DFT 
result in the original non-diagonal basis: 
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1 z 1 - iz 2 -z 3 + iz 4 
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in agreement with the original DFT transform, F m =W m „f”. 

According to various embodiments, the diagonal form of 45 
the DFT can illustrate that the DFT calculation can be per- 
formed as a vector product. This means that the DFT can be 
reduced to N multiplications of the data with the eigenvalues. 
Certain signal processing applications spend a majority of 
their time in the forward and inverse transforms part of the 50 
calculations (see FIG. 1). For signal processing applications 
where the initial and final data transformations are needed 
only once, the subsequent DFT calculations can proceed as 
simply N multiplications of the entries along the main diago- 
nal. The performance gain for a single DFT calculation can 55 
then be approximately log (N) time faster with respect to the 
FFT. For example, assuming for simplicity a power of 2 for 
the data size (N=2 W ), we can estimate a performance improve- 
ment for large n to be n log (2)=10 1 , or about an order of 
magnitude faster then the equivalent FFT calculation, and for 60 
arbitrary frequency spacing. 

Other embodiments will be apparent to those skilled in the 
art from consideration of the present specification and prac- 
tice of various embodiments disclosed herein. It is intended 65 
that the present specification and example be considered as 
exemplary only. 


What is claimed is: 

1. A method for determining aberration data for an optical 
system, the method comprising: 

providing an optical system comprising a signal processor; 

collecting intensity data generated by the optical system; 

calculating a similarity transform and an inverse of the 
similarity transform, based on the collected intensity 
data; 

calculating a discrete Fourier transform matrix in a diago- 
nal basis, based on the collected intensity data; 

pre-transforming the collected intensity data using the pro- 
cessor, by multiplying the collected intensity data by one 
of the similarity transform and the inverse of the simi- 
larity transform, to generate pre-trans formed data; 

performing an iterative transform operation on the pre- 
transformed data using the processor, by multiplying the 
pre-transformed data by the discrete Fourier transform 
matrix, a predetermined amount of times, to generate 
iteratively transformed data; 

back-transforming the iteratively transformed data using 
the processor, by multiplying the iteratively transformed 
data by one of the similarity transform and the inverse of 
the similarity transform, to form aberration data, 
wherein the iteratively transformed data is multiplied by 
the similarity transform when the intensity data is mul- 
tiplied by the inverse of the similarity transform in the 
pre-transforming step, and the iteratively transformed 
data is multiplied by the inverse of the similarity trans- 
form when the intensity data is multiplied by the simi- 
larity transform in the pre-transforming step; and 

outputting the aberration data to a user. 

2. The method of claim 1, wherein the method further 
comprises generating correction data based upon the aberra- 
tion data, and adjusting the optical system using the correc- 
tion data. 

3. The method of claim 1, wherein the outputting com- 
prises displaying the aberration data to the user. 

4. The method of claim 1, wherein the optical system 
further comprises a control unit and a telescope, the telescope 
comprises one or more optical components, at least one of the 
one or more optical components is configured to move, and 
the control unit is configured to control the movement of at 
least one of the optical components. 

5. The method of claim 4, wherein the one or more optical 
components comprise at least one mirror, lens, servo mecha- 
nism, prism, beam splitter, photomultiplier tube, detector, or 
a combination thereof. 

6. The method of claim 4, wherein the signal processor is 
operably linked to the control unit. 

7. The method of claim 1, wherein the optical system 
further comprises a space telescope. 

8. The method of claim 1, wherein the aberration data 
pertains to at least one of a piston aberration, a tip aberration, 
a tilt aberration, a defocus aberration, an astigmatism aberra- 
tion, a coma aberration, a spherical aberration, and a trefoil 
aberration. 

9. The method of claim 1, wherein the iterative transform 
operation comprises iteratively performing a phase retrieval 
operation on the pre-transformed data. 

10. The method of claim 9, wherein the method further 
comprises generating phase correction data, and the method 
further comprises adjusting the optical system using the 
phase correction data. 

11. The method of claim 1, further comprising launching 
the optical system into space prior to the collecting intensity 
data generated by the optical system. 
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12. An optical system comprising a signal processor and a 
control unit, the control unit being operably linked to the 
signal processor, wherein the signal processor is configured 
to perform a method comprising: 

collecting intensity data generated by the optical system; 5 
calculating a similarity transform and an inverse of the 
similarity transform, based on the collected intensity 
data; 

calculating a discrete Fourier transform matrix in a diago- 
nal basis, based on the collected intensity data; 
pre-transforming the collected intensity data using the pro- 
cessor, by multiplying the collected intensity data by one 
of the similarity transform and the inverse of the simi- 
larity transform, to generate pre-transformed data; 
performing an iterative transform operation on the pre- 
transformed data using the processor, by multiplying the 
pre-transformed data by the discrete Fourier transform 
matrix, a predetermined amount of times, to generate 
iteratively transformed data; 

back-transforming the iteratively transformed data using 
the processor, by multiplying the iteratively transformed 
data by one of the similarity transform and the inverse of 
the similarity transform, to form aberration data, 
wherein the iteratively transformed data is multiplied by 
the similarity transform when the intensity data is mul- 
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tiplied by the inverse of the similarity transform in the 
pre-transforming step, and the iteratively transformed 
data is multiplied by the inverse of the similarity trans- 
form when the intensity data is multiplied by the simi- 
larity transform in the pre-transforming step; and 
outputting the aberration data to a user. 

13. The optical system of claim 12, wherein the signal 
processor is further configured to generate correction data 
based upon the aberration data. 

14. The optical system of claim 13, wherein the control unit 
is configured to adjust the optical system based upon the 
correction data. 

15. The optical system of claim 12, further comprising a 
telescope, wherein the telescope comprises one or more opti- 

15 cal components. 

16. The optical system of claim 15, wherein the one or 
more optical components comprise at least one of a mirror, 
lens, servo mechanism, prism, beam splitter, photomultiplier 
tube, detector, or and combination thereof. 

17. The optical system of claim 12, wherein the optical 
system comprises a focusing system. 

18. The optical system of claim 12, wherein the optical 
system comprises a beam-generating device and a sighting 
system for the beam generating device. 



